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ABSTRACT 

We discuss an extended set of Tree+SPH simulations of the formation of clusters of galaxies, 
with the goal of investigating the interplay between numerical resolution effects and star- 
formation/feedback processes. Our simulations were all carried out in a concordance ACDM 
cosmology and include radiative cooling, star formation, and energy feedback from galactic 
winds. The simulated clusters span the mass range M^ir — (0.1 — 2.3) x 10^^ with 

mass resolution varying by several decades. At the highest achieved resolution, the mass of 
gas particles is rugas — 1-5 x 10^ which allows us to resolve the virial region of a 

Virgo-like cluster with more than 2 million gas particles and with at least as many dark-matter 
(DM) particles. Our resolution study confirms that, in the absence of an efficient feedback 
mechanism, runaway cooling leads to about 35 per cent of baryons in clusters to be locked 
up in long lived stars at our highest resolution, with no evidence of convergence. However, 
including feedback causes the fraction of cooled baryons to converge at about 15 per cent 
already at modest resolution, which is much closer to the typical values inferred from obser¬ 
vational data. Feedback also stabilizes other gas-related quantities, such as radial profiles of 
entropy, gas density and temperature, against variations due to changes in resolution. Besides 
effects of mass resolution, we also investigate the influence of the gravitational force softening 
length, and that of numerical heating of the gas induced by two-body encounters between DM 
and lighter gas particles. We also show that simulations where more DM than gas particles 
are used, such that rUgas — pidm, show a significantly enhanced efficiency of star formation 
at 3, but they accurately reproduce at z = 0 the fraction of cooled gas and the thermo¬ 
dynamical properties of the intra-cluster gas. Our results are important for establishing and 
delineating the regime of numerical reliability of the present generation of hydrodynamical 
simulations of galaxy clusters. 
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1 INTRODUCTION 

Within the hierarchy of cosmic structures, clusters of galaxies mark 
an interesting transition between the large-scale regime, where the 
dynamics is dominated by gravity, and the small-scale regime, 
where complex astrophysical processes, such as gas cooling, star 
formation and energetic feedback processes, control the formation 


and evolution of galaxies (e.g. iRosati etalJl2Q0irVoiJl2005l for 
recent reviews). If gravity were the only player at the scales rele¬ 
vant for galaxy clusters, then the overall thermal content of baryons 
would be completely determined by the processes of adiabatic com¬ 
pression and shock heatin g within the dark-matter (DM) domi - 
nated potential wells (e.g. lKaiseilll98fl : iTozzi & NormanlEoOlh . 
Since gravity has no characteristic scales, this scenario predicts 
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that galaxy clusters and groups of different mass should appear as 
scaled versions of each other, with only weak residual trends due to 
a variation of the halo concentration with mass. However, a num¬ 
ber of observational facts highlight that the evolution of cosmic 
baryons within clusters must be influenced by physical processes 
other than gravity. First, scaling relations between X-ray observ¬ 
ables demonstrate that groups and poor clusters have a lower con¬ 
tent of diffuse hot gas than rich clusters, and that they lie on a com- 


parativelv higher adiabat (e.g.lMarkevitchlll998l:lArnaud & Evrarc 

199? 

: lEttori et alJl2002l; IPonman et alJ 

l2003l lOsmond & Ponman 

2004 

). Second, the bulk of the gas in 

central cluster regions lies 


at a temperature which is never observed to fall below 1/3-1/2 of 
the virial temperature, despite the fact that th e cooling time of this 
gas is much shorter than the Hubble time (e.g .lPeterso n et all200ll: 

iMolendi & Pizzolatol200lLlBdhringer et alJ2002l) . Third, the frac- 

tion of baryons in the stellar phase with in clusters is consis tently 
rather small, around ~ 10 per cent (e.g. iBalogh et alJl200lll . with 
a tendency for poor clusters and groups to have a some what higher 
amount of stars te.g. lLin rt alfcOO-lli lR ines et a ll2Q04 . 

These observational evidences indicate the existence of a del¬ 
icate interplay between the physical processes which determine 
the evolution of the inter-galactic (IGM) and intra-cluster (ICM) 
media. Including hydrodynamical effects beyond ordinary gas dy¬ 
namics, such as gas cooling, star formation and associated feed¬ 
back processes, in simulation models is however met with sub¬ 
stantial difficulty. While analytical and semi-analytical calcula- 


tions 

(e.g. 

can provide ver 

f useful guidelines for the expected effects 

ITozz 

& NormanI 

bOOlllBower et alJl200lUBabuI et alJl2002l: 

IVoit et alJ 

20031 iLani 

et alJl2005^. we nevertheless rely on such 


simulation models, because the non-linearity of the physics and the 
geometrical complexity of cluster dynamics can usually only be 
captured in full by direct simulation. This makes it extremely im¬ 
portant to obtain an in-depth understanding of the robustness of the 
numerical methods employed, of the assumptions behind them, and 
of the sensitivity of numerical predictions for observable quantities 
on the adopted numerical parameters and approximations. 

One important validation test is to compare different simula¬ 
tion codes and different schemes for solving hydrodynamics with 
each other. For a single code, additional control tests that scrutinize 
the robustness of the predictions against resolution and different 
implementations of sub-resolution processes are highly desirable 
as well. An example for the importance of such tests is represented 
by radiative cooling, which is well known to exhibit a runaway 
character: cooling leads to an increase of the gas density, which, 
in turn, increases the cooling efficiency. Consistently, hydrody¬ 
namical cosmological simulations have shown that the fraction 
of gas that cools down and becomes available for star formation 


is much larger than indicated by observations (e.g. iKatz & White 

1993 

Suginohara & Ostrikeill99i4IPearce et alJ200fil:lLewis et al 

200C 

Yoshida et alJl2002n. a tendency that increases with resolu- 

tion ("e.g. iBalogh et alJ 

l200ll iTornatore et al.ll200.3|i. Solving this 


problem apparently requires the introduction of suitable feedback 
physics which reduces the fraction of cooled gas and stabilizes its 
value against numerical resolution. 

So far, detailed comparisons between different cosmologi¬ 
cal hydrodynamical codes h ave been restricted to the case of 
non-radiative simu lations ie.g. lKang et alll994HFrenk et alll999l : 
lOTShea et al]|2005l) . While these comparisons have shown a rea¬ 
sonable level of agreement, some sizeable differences have been 
found in the profiles of thermodynamic quantities (e.g., entropy) 
when simulations of galaxy clusters performed with E ulerian and 
Lagrangian codes were compared iFrenk et'^ll999|l . A limited 


number of studies have been presented so far which were aimed 
at discussing the effect of resolution and different implementa¬ 
tions of cooling/star-form ation physics within a single code. For 
instance, lKavetal.l f 2002h compared different implementations of 
star formation. They concluded that different prescriptions provide 
broadly consistent results, although the amount of cooled gas is al¬ 
ways much higher than observed. They also verified that this over¬ 
cooling can be partly ameliorated by resorting to either kinetic or 
thermal feedback. 

Besides checking the influence of different parameterizations 
of physical effects, it is crucial to have purely numerical effects well 
under control, such as those originating from mass resolution and 
force softening. Detailed studies have been presented about the de¬ 
termination o f an optimum softening for pur ely collisionless sim¬ 
ulations (e.g. IPower et ^ 120031: Eh^l200.4 . usually based on a 
compromise between the desire of obtaining an accurate estimate 
of the accelerations together with high spatial resolution, and at the 
same time, the need to suppress two-body relaxation and the for- 
mati on of bound particle pairs if too small a softening is chosen 
te.g. lThomas & Couchmarill992h . 

The situation is more complicated for hydrodynamical simu¬ 
lations, where energy can be spuriously transferred from the col¬ 
lisionless to the colli sional component, thereby affecting the evo¬ 
lution of the gas. As ISteinmetz & Whitel il994 have shown, this 
can happen when the gas particles are substantially lighter than the 
dark matter particles, such that the former receive a systematic en¬ 
ergy transfer in two-body encounters, leading to artificial heating 
of the gas. As a result, a lower-limit exists for the required mass 
resolution in order to provide a correct description of gas cooling 
within halos of a given mass. 


In a series of papers le.g.. iBorgani et al.l 

l2004 iDoIag et alJ 

2004 iMurante et a 

lJl2004 lEttori et al.ll2004 iDiaferio et al.ll200.5l 

Cheng et alJ200.5l 

Ettori et alJ200.5n we have presented results on 


the properties of galaxy clusters, extracted from a la rge-scale sim¬ 
ulation and performe d with the GADGET-2 code iSoringel et alJ 
l200lt I^ ringelll2005l), including the star formation and feedback 
model of lsniringd & Hemauisa i200.3alljl . This model is formu¬ 
lated as a sub-resolution model to account for the multiphase na¬ 
ture of the interstellar medium (ISM), and contains a phenomeno¬ 
logical model of energy feedback from galactic winds triggered by 
supernova (SN) explosions (see Section |3. The present paper is 
specifically aimed at discussing the numerical robustness of the re¬ 
sults we obtained, by analysing an extended set of re-simulations 
of galaxy clusters, spanning a fairly large range both in cluster mass 
and numerical resolution. More specifically, we will primarily ad¬ 
dress the following two questions: (a) How does numerical reso¬ 
lution affect the properties of the diffuse baryons and the distribu¬ 
tion of star formation within clusters? (b) How large is the impact 
of artificial heating and how can numerical parameters (i.e., force 
softening and mass-ratio between gas and DM particles) be chosen 
to minimize its effect? 

The plan of the paper is as follows. In Section 2, we provide 
the details of the simulated clusters. We analyze in Section 3 the 
effects of changing resolution over a fairly large range, up to a 
factor 45 in particle masses. In this section we will also discuss 
how the adopted feedback from galactic winds affects the resolu¬ 
tion dependence of measured properties for the simulated clusters. 
In Section 4, we discuss the role of different sources of numerical 
heating. Finally, we will summarize our results and draw our main 
conclusions in Section 5. 
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Cluster run 

Mvir 

R-vir 

Tmw 

Set#l 

CLl 

13.7 

2.3 

5.7 

CL2 

2.9 

1.4 

2.3 

CL3 

2.4 

1.3 

2.1 

CL4 

1.6 

1.1 

1.7 

Set #2 

CL5 

14.9 

2.4 

7.0 

CL6 

1.1 

1.0 

1.4 


Table 1. Basic properties of the 6 clusters belonging to the two sets of sim¬ 
ulations. Column 1: cluster name; Column 2: virial mass, defined as the 
total mass contained within the virial radius (units of 10^‘^h“^MQ); Col¬ 
umn 3: virial radius, defined as the radius encompassing an average density 
equal to the virial density predicted for the assumed cosmology (see also 
text; units of /i~^Mpc); Column 4: mass-weighted temperature computed 
within /?vir (units of keV). All the values reported refer to the reference 
runs of both simulation sets (see text). 

2 THE SIMULATIONS 

Our simulations were carried out with GADGET-2 jSDringell2003). 
an impro ved version of the p arallel Tree-SPH simulation code 
GADGET iSpringel et all20 ()In. It uses an entro py-conserving for¬ 
mulation of SPH Isnrin^ & Hemnuislll20o3i . and includes ra¬ 
diative c ooling, heating by a uni form redshift-evolving UV back¬ 
ground jHaardt & MadaJll996t . and a treatment of star forma¬ 
tion and feedback processes. The prescription of star formation is 
based on a sub-resolution model to account for the multi-phase 
nature of the interstellar medium (ISM), w here the cold phase of 
the ISM is the reservoir of star forming gas iSoringel & Hemauisj 
l2003aL SH03 hereafter). As for the feedback, SH03 followed a 
phenomenological scheme to include the effect of galactic winds, 
whose velocity, Vw, scales with the fraction rj of the SN-II feedback 
energy that contributes to the winds, as oc (see eq.[28] in 
SH03). The total energy provided by SN-II is computed b y assum¬ 
ing tha t they originate from stars with mass > 8 Mq for a lSalpeteJ 
^19■5.5^ initial mass function (IMF), with each SN releasing 10®^ 
ergs. As discussed in the following, we will assume 77 = 0.5 and 
1, yielding Vw — 340 and 480 km s“^, respectively, while we will 
also explore the effect of switching off galactic winds altogether. 

We consider two sets of clusters, which have been selected 
from different parent cosmological boxes. Initial conditions for 
both sets have been ge nerated using t he Zo omed Initial Condi¬ 
tion (ZIC) technique by iTormen et alJ il997h . This technique in¬ 
creases the mass resolution in a suitably chosen high-resolution 
Lagrangian region surrounding the structure to be re-simulated. It 
then adds additional initial displa cements, assigned according t o 
the Zeldovich approximation fe.g. IShandarin & Zeldovict]ll989ll . 
from the newly sampled high-frequency modes which were not 
assigned in the low-resolution parent simulation. Furthermore,the 
mass resolution is progressively degraded in more distant regions, 
so as to save computational resources while still correctly describ¬ 
ing the large-scale tidal field of the cosmological environment. 

Once initial positions and velocities are assigned, a DM-only 
run is performed to check for any contamination of the surround¬ 
ings of the cluster virial region by heavy particles that may have 
moved in from the low- to the high-resolution region. If required, 
the shape of the Lagrangian high-resolution region is optimized by 
trial and error until any such contamination around the halo of inter¬ 
est is prevented. With a typical number of 3-5 trials, we end up with 


Run name 

mDM 

JTlgas 

epi 


Set #1 LR 

46.4 

6.93 

7.5 

2 

Set #1 MR 

15.5 

2.31 

5.2 

2 

Set #1 HR 

4.6 

0.69 

3.5 

2 

Set #1 VR 

1.0 

0.15 

2.1 

2 

Set #2 

11.3 

1.69 

5.0 

5 


Table 2. Parameters defining the two sets of runs (see also text). Column 
1: name of the run; Column 2: mass of the high-resolution DM particles 
(in units of 10 ®/i“^Mq); Column 3: mass of the gas particles (in units of 
lO®/i“^M 0 ); Column 4: Plummer-equivalent gravitational force soften¬ 
ing at 2 = 0 in the high-resolution region (units of fi~^kpc). Column 5: 
redshift of transition from physical to comoving softening. 


initial conditions which produce a cluster that, at 2 = 0 , is free of 
contaminants out to (4-6) Rvii- Once initial conditions are created, 
we split particles in the high-resolution region into a DM and a gas 
component, whose mass ratio is set to reproduce the assumed cos¬ 
mic baryon fraction. Instead of placing them on top of each other, 
we displace gas and DM particles such that the centre of mass of 
each parent particle is preserved and the final gas and dark matter 
particle distributions are interleaved by one mean particle spacing. 


2.1 The Set #1 

This set includes four clusters, resimulated at different resolutions, 
with virial mass in the range Mvir = (1.6-13) x 10 ^"*/i“^Mq (CLl 
to CL4 in Table0. These clusters have been extracte d from the cos¬ 
molog ical hydrodynamical simulation presented bv ISorgani et alJ 
i2004ll . The simulation followed 480® DM particles and an ini¬ 
tially equal number of gas particles, within a box of 192 /i“®Mpc 
on a side, for a flat ACDM model with fl™ = 0.3, h = 0.7, 
(78 = 0.8 and fib = 0.04. With these choices, the masses of 
the DM and gas particles are ttidm — 4.6 x 10® /i~®M 0 and 
rrigas — 6.9 x lO®/i“®M 0 , respectively. The force accuracy is 
set by epi = 7.5/i“®kpc for the Plummer-equivalent softening 
parameter, fixed in physical units from 2 = 0 to 2 = 2 , and kept 
fixed in comoving units at higher redshifts. 

Since initial conditions of the parent simulations have been 
generated on a grid, a grid is also used to assign initial displace¬ 
ments for these resimulations. Initial conditions for each cluster 
are generated at four different mass resolutions, corresponding to 
the basic resolution of the parent box (low resolution, LR), and to 
3 times (medium resolution, MR), 10 times (high resolution, HR) 
and 45 times (very high resolution, VR) smaller particle masses. 
The VR run is not carried out for the CLl cluster, whose large 
mass would result in a too large computational cost. The gravi¬ 
tational softening in the high-resolution regions is rescaled with 
the mass of the particles according to epi oc where the LR 

runs were set to have the same softening as used in the parent sim¬ 
ulation. In Table 1^ we list the masses of the DM and gas parti¬ 
cles, as well as the force softening for the different resolutions. 
At the highest achieved resolution, each cluster is resolved with 
at least 1,5 million DM particles within the virial radius, with CL2 
reaching 2,2 million particles. The reference runs for the clusters 
of Set #1 have been performed by assuming that 100 per cent of 
the energy provided by SNe is carried by winds. This gives a wind 
speed of Vw ~ 480kms“®. 
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Run name 

Description 

Applied to: 

NW 

No winds: winds switched off 

CL2-LR,HR; CL4-LR,MR,HR,VR 

sw 

Strong winds: wind velocity for the Set #2 at the reference value of the Set #1 

CL5, CL6 

108 

Eight times heavier gas particles 

CL4-VR 

SOS 

SW runs for clusters of the Set #2 with ug = 0.8 

CL5, CL6 


Table 3. Description of the name extensions used for the identification of the runs. 


2.2 The Set #2 

This set includes 2 clusters having mass ~ 10^® and 

~ 10^“* h“^M0 (CL5 and CL6 in TableQ respectively). They have 
been selected from a larger sample of 20 clusters with masses in the 
range 5 x 10^® — 2.3 x 10^® which have been identified in 

simulations of 9 Lagrangian regions (Dolag et al., in preparation). 
These systems were extracted from a DM-only simulation with a 
box size of 479 /i“^Mpc of a flat ACDM model with Om = 0.3, 
h = 0.7, (Tg = 0.9 and fib = 0.04 iYoshida et alJl20()Ill . Differ¬ 
ently fr om Set #1, th e initial displacements were generated using a 
‘glass’ JWhita 199611 for the Lagrangian particle distribution. Only 
one mass resolution was simulated in this case, corresponding to 
muM = 1-13 X 10® h“^M0 and mgas = 1.7 x 10® /i“^M0 for 
dark matter and gas within the high-resolution region, respectively. 
As such, this mass resolution is about 4 times better than the LR 
runs, and therefore lies intermediate between the MR and the HR 
runs of Set #1. The softening length is set to epi = 5.0/i“^kpc, 
fixed in physical units below z — 5, while it is kept fixed in co¬ 
moving units at higher redshift. The reference runs for clusters of 
Set #2 assume that 50 per cent of the energy provided by SN (i.e. 
half of that used for the Set #1) is carried by winds. This gives 
a wind speed of Vw — 340kms“^. Unless otherwise stated, our 
analysis presented in this paper is restricted to the CL5 and CL6 
simulations out of the 20 clusters of the full Set #2. We note that 
for all of our simulations, we set the smallest allowed value for the 
SPH smoothing length to epi/4. 

In summary, our two sets of simulated clusters differ in the 
following aspects; 

(a) The value of erg changes from 0.8 to 0.9 when going from 
Set #1 to Set #2. To explicitly check for the effect of changing the 
normalization erg, we will also analyse additional simulations of the 
CL5 and CL6 clusters of Set #2 using erg = 0.8. 

(b) The initial Lagrangian positions of DM and gas particles were 
taken to be a grid for Set #1 and a glass for Set #2. We will also 
present results for CL5 of Set #2 using grid initial conditions, to 
check explicitly for effects due to a grid vs. a glass setup. 

(c) The velocity of the galactic winds is n™ ~ 480kms“^ and 
340 km for Set #1 and Set #2, respectively. To address this im¬ 
portant difference, we have rerun the CL5 and CL6 clusters from 
Set #2 with erg = 0.8 and with the stronger feedback, to make 
them fully comparable with the clusters of Set #1. 

(d) The transition from physical to comoving softening takes 
place at z = 5 for Set #2 and at 2 = 2 for Set #1. To test the 
effect of a change of this transition redshift we have rerun CL6 of 
Set #2 with strong winds, also using the same choice of softening 
as for Set #1. 

We will use the convention that a certain simulation of each 
cluster is labeled by the name of the cluster itself, followed, when 
required, by a label which specifies the resolution used. In this way, 
CL I-MR will indicate the reference run of medium-resolution of 


the CLl cluster from Set#l, while CL5 will designate the reference 
run of this cluster from Set #2. For the latter, we do not specify 
the resolution, since clusters from Set #2 are simulated at only 
one resolution. In addition, a further extension of the name of each 
run will be provided whenever the run differs from the reference 
run of the set it belongs to. The list of such extensions and their 
description is given in Tabled For instance, CL4-HR indicates the 
high-resolution run of the CL4 cluster, using the standard setup 
of Set #1, while CL4-HR-NW stands for the same simulation, but 
neglecting the effect of galactic winds. 

We show in FigureQthe gas density maps of the six simulated 
clusters within boxes each having a size of 4i7vir- The four clusters 
from the Set #1 are shown in their highest resolution version. 


In the first step of our simulation analysis we determine suitable 
cluster centers. These are defined as the position of the most bound 
particle among those grouped together by a FOF algorithm with 
linking length b = 0.15 (in units of the mean interparticle sepa¬ 
ration in the high-resolution region). Once the center is identified, 
we apply a spherical overdensity algorithm to determine the virial 
radius, Rvii- We here defined this as the radius that encompasses an 
average density equal to the virial density for the adopted cosmo¬ 
logical model, = Ac{z)pciz) (Pc{z) = [H(z)/Hof Pc,o 

is the critical density at red shift 2), where th e overdensity Ac (2) 
is computed as described bv lEke eta i] jl996h . with Ac(0) ~ 100 
for the assumed cosmology. The virial mass, Mvir, is simply the 
mass contained within the virial radius. In TableQ we provide the 
typical values of Rvii, Mvir and mass-weighted temperature, Tmw, 
for the six simulated clusters. Profiles of gas-related quantities are 
computed within 200 equispaced linear radial bins, out to IR^ii, 
starting from a minimum radius w hich contains 100 gas particles. 

As shown bv iBorgani et al ] <2002h . numerically stable results can 
be expected for this choice in non-radiative simulations of galaxy 
clusters. 

We ide ntify galaxie s in our simulations by applying the SKID 
algorithm^ iStadell200lh to the distribution of star particles. In the 
following, we provide a short description of our algorithm, while a 
more detailed discussion and presentation of tests is provided else¬ 
where (Murante et al. 2005, in preparation). Briefly, the SKID al¬ 
gorithm works as follows: 

• An overall density field is computed by using the distribution 
of all the particle species, i.e. DM, gas and star particles. The den¬ 
sity is estimated with a SPH spline-kernel, using a given number 
Nsrn of neighbour particles. 

• The star particles are moved along the gradient of the density 
field in steps of r/2. When a particle begins to oscillate inside a 
sphere of radius t/2, it is stopped. In this way, r can be interpreted 

^ See http://www-hpcc.astro.Washington.edu/tools/skid.html 


2.3 Basics of the simulation analysis 
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Figure 1. Maps of the gas density for the six simulated clusters, in a region encompassing two virial radii around each cluster centre. Top-left to bottom-right 
panels are for the CLl to CL6 clusters. For the clusters belonging to Set #1 (CLl to CL4), we show maps for the highest resolution runs (HR for CLl and VR 
for CL2 to CL4). In all cases, the maps refer to the runs with strong winds (i.e. Vw — 480 km s~^). The small bright knots mark the “galaxies”, i.e. the places 
where high density gas is undergoing cooling and star formation. 
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Figure 2. The fraction of cooled baryons /c (left panel), and the number of galaxies within Vir, as functions of the mass of the gas particle, for the clusters 
of Set #1 at different resolution. Filled symbols are for the “reference” runs. The open circles are for the runs of the CL4 cluster with wind feedback turned 
off. The asterisk is for the CL4 run at very high resolution using 8 times heavier gas particles (CL4-lo8), so that the gas particle mass is similar to that of the 
DM particles in the high-resolution region. 


as the typical size of the smallest resolved structure in the distribu¬ 
tion of the star particles. 

• When all particles have been moved, they are grouped using 
a friends-of-friends (FOF) algorithm applied to the moved particle 
positions. The linking length is again t/2. 

• The binding energy of each group identified in this way is 
computed by accounting for all the particles inside a sphere cen¬ 
tered on the center of mass of the group and having radius 2r (for 
the moved particles, their initial positions are used in the computa¬ 
tion of the potential). This binding energy is then used to remove 
from the group all the star particles which are recognized as un¬ 
bound. Finally, we retain such a SKID-group of stars as a galaxy if 
it contains at least 32 particles after the removal of unbound stars. 

The resulting list of objects identified by SKID depends on 
the choice of two parameters, namely r and Nsm- After many ex¬ 
periments, and resorting to visual inspection in all cases, we found 
that a complete detection of bound stellar objects requires use of 
a set of different values of Nsm- We used Nsm = 16, 32, 64, and 
define a “galaxy” to be the set of star particles which belong to a 
SKID group with any one of the above Nsm values. If a star parti¬ 
cle belongs to a SKID group for one value of Nsm and to another 
group for a different Nsm-, then the groups are “joined” and are 
considered as forming a single galaxy. All star particles not linked 
to any galaxy are considered to be part of a diffuse stellar compo¬ 
nent in the cluster iMurante et alJ2004 . As for T, since it roughly 
corresponds to the size of the smallest resolved structure, we adopt 
r = 2.8epi, which is the scale where the softened force becomes 
equal to the Newtonian force. 


3 THE EFFECT OF MASS RESOLUTION 

Checking the stability of simulation results against numerical res¬ 
olution is always of paramount importance to assess their robust¬ 
ness. For hydrodynamical simulations of clusters that only account 


for non-radiative physics, it has been shown that numerically reli¬ 
able estimates of global cluster properties, such as temperature and 
X-ray luminosity, are obtained when each halo is resolved with a 
few tens of thousands of gas particles within the virial region, a 
resolution which is also enough to result in converged radial pro¬ 
files for the gas density, tempe rature and entropy down to a few 
perce nt of the virial radius (e.g., iNavarro et alll99^lBorgani et alJ 
I2OO2I and references therein). However, the situation is consider¬ 
ably less clear when additional physical processes are introduced. 
For example, the efficiency of radiative cooling in simulations is 
known to sensitively depend on the adopted mass resolution. At the 
same time, the effects of feedback from galactic winds may also 
depend on numerical resolution, and it is unclear whether this reso¬ 
lution dependence favourably counteracts the increasing efficiency 
of cooling. We generated the initial conditions for Set #1 at four 
different resolutions precisely for the purpose of checking in detail 
the resolution-dependence of the non-trivial physical processes in¬ 
cluded in our simulations. 


3.1 Star formation 

In Figure we show the dependence of the fraction of cooled 
baryons and of the number of identified galaxies within iivir on 
the mass of the gas particles. In the absence of feedback by galac¬ 
tic winds, the CL4 runs develop the typical runaway of cooling as 
a function of the mass resolution. The fraction of cooled baryons 
steadily increase from 21 per cent at the lowest resolution to 34 per 
cent at the highest achieved resolution, with no indication of con¬ 
vergence. A part of the origin of this runaway is illustrated in the 
left panel of Figure|3 where we show the star formation history for 
the CL4 runs at the different resolutions, in the absence of winds. 
Increasing the resolution leads to an earlier onset of star formation 
within ever smaller first-collapsing halos, and at low redshift, cool¬ 
ing and star formation appear to never enter in any self-regulated 
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Figure 3. The star formation rate as a function of the mass resolution for the CL4 cluster, both excluding (left panel) and including (right panel) the effect 
of galactic winds. In both panels, solid, short-dashed, long-dashed, and dotted curves correspond to the VR, HR, MR and LR runs, respectively. In the right 
panel, the heavy solid curve is for the VR run, where gas particles have nearly the same mass as the DM particles (lo8 run in Table|^. Note the different scales 
on the vertical axes in the two panels. 


regime. As a result, the star formation efficiency progressively in¬ 
creases with resolution at all epochs. 

On the other hand, including feedback by galactic winds is 
quite effective in regulating the process of star formation. In this 
case, there is no evidence for a systematic increase of the fraction 
of cooled baryons with increasing resolution. In fact, the three poor 
clusters have rather similar values of ~ 15 per cent, and the 
massive cluster has a value of fc — 12 per cent. This result demon¬ 
strates that a converged estimate of the star fraction is obtained al¬ 
ready at modest resolution with our feedback scheme. Th is also 
extends the earlier results bv ISpringel & Hernouisll i2003lj) to the 
scale of clusters. Intriguingly, we note that all the runs at very high 
resolution (VR) even show a small but systematic decrease of the 
cold gas fraction. 

Including feedback from winds has the twofold effect of sup¬ 
pressing the resolution-dependence of the cumulative efficiency 
and to make the SFR history almost independent of resolution at 
zf} 2-3 (see the right panel of Fig.|3- The enhanced star formation 
activity at high-redshift of the high resolution runs produces more 
efficient gas pre-heating, which has the effect of inhibiting star for¬ 
mation at later times. This presumably explains that the star forma¬ 
tion rate of the VR run at 2 is even below those of the lower- 
resolution runs. Assuming that the smallest resolved halos where 
cooling can take place contain ~100 DM particles, we find their 
escape velocity to be of order lOOkms”^ for the low-resolution 
(LR) runs, and a factor of about three smaller for the VR runs of 
Set #1. Therefore, the wind velocity is always larger than the es¬ 
cape valocity of the smallest halos which are first resolved at high 
redshift, thereby implying that our feedback scheme is efficient in 
preventing star formation already in the first generation of resolved 
galaxies. 

An alternative explanation for the decrease of the cold frac¬ 


tion in the VR runs could be that the suppression of star forma¬ 
tion in the highest resolution run is related to some numerical ef¬ 
fect. For instance, one may argue that an improved resolution pro¬ 
vides a more accurate description of the gas behaviour at cooling 
interface. Indeed, a coarse desc ription of this interf ace is expected 
to cause spurious gas cooling iPearce et alJIlQQ^^l . a feature that 
shoul d however be weak in the en tropy conserving formulation of 
SPH iSoringel & HemauisJ20o3) implemented in our code. How¬ 
ever, if this was really the case, one would expect the same effect to 
appear also in runs without galactic winds (NW runs). But the left 
panel of Fig.l^clearly demonstrates that this is not the case. Hence, 
the stable behaviour of star formation with increasing resolution is 
more likely related to the inhibiting effect of more efficient high-2 
feedback onto later generations of galaxies. 

The number of identified galaxies within each cluster in¬ 
creases with better resolution and with the cluster mass, as ex¬ 
pected. As shown in the right panel of Figure the number of 
identified bona-fide galaxies grows by more than one order of mag¬ 
nitude when passing from the LR to the VR runs. In line with the 
behaviour of the stellar fraction, the number of galaxies for the CL4 
cluster increases by a factor of about three when the wind feedback 
is switched off, almost independent of the resolution. We defer a de¬ 
tailed analysis of the stellar mass function of the identified galaxies 
and of the diffuse stellar component as a function of resolution to 
forthcoming work (Murante et al., in preparation). 

In Figure |4| we summarize the results on the star fraction by 
plotting it as a function of the cluster virial mass, for the reference 
runs of both Set #1 (filled circles) and of all the 20 clusters of 
Set #2 (squares). Both simulation sets confirm a trend of decreas¬ 
ing star fraction as a function of the cluster mass. However, clusters 
belonging to Set #2 show a stellar fraction which is systematically 
higher than that of Set #1, the difference being larger than any pos- 
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Figure 5. The effect of changing the power spectrum normalization cjs. Left panel: the effect on the star-formation history for the CL5 and CL6 clusters 
(upper and lower pairs of curves), both simulated with the strong-wind (SW) feedback; solid and dashed curves refer to erg = 0.9 and erg = 0.8, respectively. 
Right panel: the effect on gas-related profiles for the CL5 cluster. 





Figure 4. The stellar fraction, /*, as a function of the cluster virial mass, 
for the different simulation sets. Filled circles connected by a dashed line 
are for the high-resolution version of the clusters of Set #1. Squares are for 
all the clusters of Set #2, with the two filled squares showing the results 
for simulations CL5 and CL6. Filled triangles ai‘e for CL5 and CL6, but 
using the same velocity of galactic winds, v-w ~ 480 km s ^ as for Set #1. 
The open triangles are the same simulations as the filled triangles, but using 
the power spectrum normalization of Set #1,0-8 = 0.8. Therefore, filled 
circles and open triangles come from different sets of initial conditions, 
which however have been evolved using an identical simulation set-up. 

sible object-to-object intrinsic scatter induced by the varying dy¬ 
namical histories of different clusters. However, we recall that the 
two sets of cluster simulations differ in the strength of the adopted 
feedback and in the normalization of the power spectrum. 


To investigate the effect of the different feedback, we have run 
additional simulations of CL5 and CL6 of Set #2 by increasing 
the wind speed to the same value as used for Set #1 (SW runs). 
The results for the stellar fraction are shown with filled triangles 
in Fig.|3 Although increasing the feedback efficiency produces a 
significant suppression of /,, the effect is still not large enough to 
fully account for the difference between the two simulation sets, 
thus suggesting that the residual difference is due to the different 
(78 values. 

We have explicitly verified this by repeating the CL5-SW and 
CL6-SW runs by also decreasing as to 0.8. The comparison of 
the resulting star formation histories is shown in the left panel of 
Figure|3 Reducing as results in a significant change in the timing 
of structure formation and, correspondingly, to a delay in the high- 
redshift star formation, with a suppression of its peak at 2 ~ 3. The 
resulting stellar fraction further decreases from /* = 17.5 per cent 
to 15.6 per cent for CL5, and from /* = 20.3 per cent to 18.1 per 
cent for CL6. As shown in Fig.|3 the reduction of /* connected 
to the power spectrum amplitude finally brings the values of the 
cooled fraction for the simulated clusters of Set #1 and Set #2 into 
good agreement. 

The decreasing trend of the star fraction with cluster mass, 
shown in Fig. IH is in qualitative agreem ent with observational 
results jLin_et_allj200j: iRines et alJl2004 cf. also iBalogh et alJ 
1200 ih . lLh^i^nll200^ used K-band data from the 2MASS sur¬ 
vey to trace the stellar population in clusters, and ROSAT-PSPC 
data to measure the corresponding gas mass. They found /* ~ 10 
per cent for clusters with mass of about 10^® increasing 

to ~ 15 per cent for clusters with 10^"^ These values are 

not far from those obtained from our runs with stronger feedback 
and as = 0.8. Also, it is worth pointing ou t that the observa tional 
census of the cluster stellar population by iLin et all i200.3h does 
not include diffuse stars, whose cont ribution may not be negligible 
both for real (e.g., J^niaboldi||2004 , for a review) and simulated 
iMurante et al.l l2004 Iwillmm^^ll 120041 : ISommer-Larsen et"^ 
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Figure 6. The effect of resolution on the profiles of gas density, temperature and entropy (from top to bottom in each panel). Each panel shows the profiles for 
the different resolutions at which each cluster of Set #1 has been simulated. For the CL4 cluster, we also show the results for the series of runs where the effect 
of winds has been excluded (NW). 
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Figure 7. The distribution of gas particles for the very high-resolution (VR) version of the CL4 cluster, lying within O.lRvir around the centre, in the rie-T 
plane (rie and T: electron number density and temperature of each gas particle, in units of cm“® and keV, respectively). Left and right panels are for the run 
without and with the effect of galactic winds, respectively. In each panel, the dashed curves indicate levels of constant entropy {S = 1000, 300, 100, and 30 
keV cm^, from upper to lower curves). 


|2005^ clusters. Note that our estimates of /* include all stars found 
inside the simulated clusters. 

How do our results compare with previous simulations? 
An overproduction of stars in simulations including radiative 
physics has been a well established result for sever al years (e.g., 
Katz & Whitl 1 199.31: Suginohara & Ostrikeil 1 19981 : iLewis el's] 
20001: iDave et alJl2002) . This problem is generally thought to be 
solved by a physical feedback process, which prevents excessive 
gas cooling with a continuous energy supply. However, realisti¬ 
cally mo delling such a h e ating process is a difficult problem. For 
instance, iTornatore et alJ ^2003^ verified that using thermal feed¬ 
back instead of the kinetic one applied here, the fraction of cooled 
baryons is of order 30 per cent, with no evidence for con vergence 
with better numerical resolution. iMuanwong et a i U2Q0g) claimed 
that overcooling within clusters can be e ventually avoided by suit¬ 
ably pre-heating gas at high redshift. iKav et all i2004l imple¬ 
mented a scheme of thermal feedback which increases the entropy 
of gas particles which are just undergoing cooling. Although this 
feedback is tuned to reproduce several X-ray observational prop¬ 
erties of clusters, it still provides a too large star fraction, of about 
25 per cent . Whil e all these results are based on SPH codes, 
iKravtsov et alJ used an adaptive Eulerian code, which also 

includes cooling, star formation and thermal feedback. They find 
rather large stellar fractions as well, with values of about 30-40 per 
cent within the cluster virial radius, which is consistent with the 
SPH simulations. 

3.2 Thermodynamics of the hot gas 

In Figure|^ we show the radial profiles of gas density, temperature, 
and entropy, for the the four clusters the Set #1, at different resolu¬ 
tions. We find that the profiles are rather stable against resolution on 
scales O.liZvir. On smaller scales, where the complex physics 
of cooling, star formation and winds’ feedback plays a significant 
role, we detect significant systematic changes of the profiles with 


increasing resolution. In general, the gas density and entropy pro¬ 
files become slightly shallower, and the temperature profiles some¬ 
what steeper, as the resolution is increased. The size of this effect 
changes from object to object, depending on the different dynam¬ 
ical histories, being more apparent for the CLl cluster and rather 
small for the CL3 cluster. 

Shallower gas-density and entropy profiles can be due ei¬ 
ther to more efficient non-gravitational heating from energy feed¬ 
back or to more efficient radiative cooling. In the first case, the 
stronger heating places gas on a higher adiabat, thereby prevent- 
ing it from reaching h igh density in the central halo region (e.g. 
iTozzi & Normanl200lL and references therein). In the second case, 
a more efficient cooling turns into a more effic ient selective remova l 
of low entropy gas from the hot phase (e.g., IVoit fc Rrva^l2n0lh . 
As a result, the gas density decreases, while leaving behind the hot 
phase at its relatively higher entropy. If cooling is the process that 
governs the resolution-dependence of the profiles, we would ex¬ 
pect a more efficient star formation and a larger fraction of cooled 
gas at increasing resolution. However, this is not found for the runs 
including galactic winds (see Figs.|^and|3- On the other hand, in¬ 
creasing the resolution has the effect of increasing star formation 
efficiency and, therefore, the amount of feedback energy released 
to the gas at high redshift. Indeed, galactic winds are expected to be 
more efficient at high redshift, when they can more easily escape 
from shallower potential wells. Better resolution hence increases 
the gas heating associated with high-redshift star formation. 

The inner slope of the temperature profiles obtained in radia¬ 
tive simulations is k nown to be significan tly stee per than observe d 
for real clusters (e.g.. lBorgani et all2004 cf. also lKav ^ all2004h . 
This effect has been interpreted as being due to the adiabatic com¬ 
pression of gas which falls in from outer cluster regions, as a con¬ 
sequence of the reduced pressure support after removal of cooled 
gas. Our results show that improving resolution does not help to 
reconcile simulated and observed temperature profiles. In fact, this 
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Figure 8. The fraction of cooled baryons (left panel) and the number of galaxies (right panel) within the virial radius, as a function of the gravitational 
softening, for the CL2 and CL4 clusters (both in the LR and HR version). The softening is given in units of the value assumed for the reference runs (see 
Tablel^. These simulations did not include galactic winds (NW rans). 


discrepancy can become even worse with better numerical resolu¬ 
tion. 

In Fig.|^ we show profiles of the gas properties for the CL4 
runs with no winds (NW). Also for these runs, we note similar 
trends with resolution as for the runs including feedback by winds. 
However, in this case the increase of the central entropy level, and 
the corresponding decrease of the gas density, are explained by the 
higher efficiency of cooling with better resolution. Therefore, al¬ 
though the runs with and without feedback by winds show similar 
trends with resolution, the interpretation of these effects is com¬ 
pletely different in the two cases. 

Furthermore, a comparison of the CL4 runs with and without 
winds shows that, at fixed resolution, the former have higher central 
density, shallower temperature profiles and lower entropy levels. In 
order to explicitly demonstrate the effect of feedback on the gas 
thermodynamics in the central cluster regions, we show in FigureQ 
a phase diagram of the gas particles lying within the central 0. li?vir 
for the CL4 cluster, with and without feedback by galactic winds. 
Perhaps somewhat counterintuitively, the effect of feedback is that 
of lowering the typical temperature of the gas in the central cluster 
regions, and shifting it to a lower adiabat. This result emphasizes 
the role of feedback in the central cluster regions: a continuous sup¬ 
ply of energy by feedback has the effect of keeping a population of 
relatively low-entropy gas particles in the hot phase, which have a 
short enough cooling time that they would otherwise have dropped 
out of the hot phase and cooled onto the centre. The presence of 
such gas particles explains the higher gas density and the lower 
entropy level. Furthermore, the higher pressure support associated 
with the feedback energy reduces the compressional heating of in¬ 
falling gas and, as a consequence, the ICM temperature in the cen¬ 
tral cluster regions. 

Observational data on the temperature structure in cool-core 
clusters indicate that the ICM reaches there a lower limiting 
temperature whic h is about 1/2-1/4 of the overall cluster virial 
temperature (e.g.. |Peterson et all200ll:lMolendi & Pizzolatol200ll: 

iBbhringer et al J2OO2I) . The generally accepted interpretation is that 


some sort of feedback prevents gas from reaching lower tempera¬ 
ture by compensating its radiative losses, thus suppressing also the 
cooling rate. Our results on the effect of feedback in simulations 
clearly corroborates this picture. However, our implemented feed¬ 
back is still not efficient enough to suppress the cooling rate and, 
correspondingly, the compressional heating to the observed level. 
As a consequence, the temperature profiles in the central regions 
are still steeper than observed, and no cool cores are created. 

Finally, we consider the effect of changing as on the profiles 
of gas-related quantities (right panel of Fig.|3. We find that the gas 
properties are left almost unchanged, once they are rescaled to com¬ 
pensate for the difference in virial temperature. However, slightly 
shallower density and entropy profiles, and a slightly steeper tem¬ 
perature profile characterize the runs with erg = 0.9. This can be 
understood as a consequence of the enhanced cooling, which re¬ 
moves a larger fraction of gas from the hot phase in the central 
cluster regions and, correspondingly, increases the compressional 
heating of the in-flowing gas. 


4 THE EFFECT OF NUMERICAL HEATING 

4.1 Changing the gravitational softening 

The choice for the gravitational softening is a compromise be¬ 
tween the goal of resolving the smallest scales possible and the 
need to ensure that the relaxation time by two-body interactions 
is much larger than the typi cal age of the simulated structure (e.g. 
iThomas & Couchmanll992h . Using a very small softening may al¬ 
low one to resolve a larger number of small halos, but this comes 
at the price of introducing spurious numerical heating of the gas, 
which could artificially suppress cooling. Therefore, by increasing 
the softening we expect to first see a more efficient cooling, until 
a point is reached when only large halos, with their comparatively 
long cooling times, are resolved and the cooling efficiency declines 
again. This implies that a softening value should exist which max- 
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imizes the amount of cooled gas, while smaller and larger values 
lead to less cold gas as a result of spurious gas heating or lack of 
resolution, respectively. 

To check for this effect, we show in Figure|^the dependence 
of the fraction of cooled haryons, /c, and of the number of iden¬ 
tified galaxies rigai within the virial radius on the adopted gravita- 
tiobal softening length. We show results both for the low-resolution 
(LR) and the high-resolution (HR) versions of the CL2 and CL4 
clusters. In order to avoid mixing the effects of numerical heating 
with that of efficient feedback, we have performed these simula¬ 
tions by switching off feedback by galactic winds (NW). In the plot 
of Figure]^ the softening is given in units of that adopted for our 
reference simulations, as reported in Table The results clearly 
confirm our expectation: both a too large and a too small softening 
lead to a decrease of the amount of cooled gas. Quite interestingly, 
the softening (in units of the reference value) at which the cooled 
baryon fraction is maximised is the same at low and high resolu¬ 
tion. This demonstrates that the adopted scaling of the softening 

1 /3 

length with mass resolution, epi oc rngL is in fact a reasonable 
choice. 

In general, we find that the maximum value of /c is attained 
for a softening value about twice as large as our reference value. 
This implies that a two times larger softening should then be pre¬ 
ferred, although the difference in /c when using epi = 1 instead of 
2 is always rather small, especially at high mass resolution. 

As for the number of identified galaxies, it monotonically de¬ 
creases with increasing softening (right panel of Fig. |^. Using 
epi = 2 provides a ~ 30 per cent smaller number of galaxies than 
for epi = 1. However, this argument does not necessarely repre¬ 
sent a valid support for the choice of a smaller softening. Indeed, 
a spurious energy transfer from DM to gas may lead to a tighten¬ 
ing of the DM halos, thereby producing a larger number of small 
halos. Clearly, a close investigation of the optimal strategy to mini¬ 
mize the effect of numerical heating would require a more detailed 
analysis than that presented here, including a close investigation of 
the mass function of the resolved galaxies. We plan to present this 
analysis in a forthcoming paper. 

We finally note that the anomalous number of galaxies for the 
HR version of the CL2 cluster is apparently a result of a misclas- 
sification by SKID, as we checked by visual inspection. Here, the 
distribution of star particles shows that they are arranged in a few 
over-merged galaxies, most of which have pretty large tidal tails. 
The SKID algorithm sometimes splits these tails into separate grav¬ 
itationally bound structures, which are then misclassified as galax¬ 
ies. 

As a further test, we have compared the results of the CL6 run 
with strong winds, assuming z = 2 and z = 5 for the redshift of 
transition from comoving to physical softening. As shown in Fig- 
urej^the smaller high-2 softening provides an enhanced efficiency 
of star formation. In this case, the cooling efficiency within small 
resolved halos dominates over the effect of numerical heating in 
defining the pattern of star formation. 


4.2 Gas-DM particle splitting 

As we have already discussed, heating of gas particles from 
two-body encounters may spuriously alter the thermodynamic 
properties of the gas. Based on an estimate of the heating 
time from two-body enco u nters in the impuls e appro ximation 
jBinnev & Tremaind ISteinmetz & Whit j Il997li analyti¬ 

cally derived the heating time associated with such encounters. Un- 



Figure 9. Comparison of the results for the CL6 cluster using both 2 = 5 
(solid line) and 2 = 2 (dashed line) for the transition redshift from physi¬ 
cal to comoving softening. Strong winds (SW) have been included in these 
simulations. 


der the assumption of hydrostatic equilibrium, this heating time is 
inversely proportional to the mass of the DM particle, while being 
independent of the gas particle mass. As a consequence, a limiting 
mass of DM particles can be derived for two-body heating rate to 
dominate over the cooling rate within a given halo, as a function 
of the virial temperature of the halo itself. According to their re¬ 
sults, the mass resolution achieved in our runs, with the possible 
exception of the low-resolution (LR) series for clusters of Set #1, 
should be adequate to prevent spurious heating within halos hav¬ 
ing temperature of at least few keV. While this holds for already 
formed halos at low redshift, the situation may be different at high 
redshift, when gas cooling takes place within small, just resolved 
galaxy-sized halos. Also, in the first generation of halos, substan¬ 
tial residual gas motions exist which could give rise to an additional 
artificial energy transfer in two-body encounters from the heavier 
DM particles to the lighter gas particles. 

In order to reduce the amount of numerical heating from such 
energy equipartition effects one may try to increase the number 
of DM particles, thereby decreasing their mass, while keeping the 
mass resolution of the gas fixed. In order to check for this effect, 
we resimulated the CL4 halo assuming very high resolution (VR) 
for the DM particles. However, in this case we generated initial dis¬ 
placements for an eight times smaller number of gas particles (lo8 
runs in Table |3. In this way, gas and DM particles have a com¬ 
parable mass: mgas/moM — 1.06 instead of ~ 0.15. In terms 
of gas-mass resolution, this simulation is intermediate between the 
MR and HR runs, but with a more accurate treatment of gravity, 
thanks to the larger number of DM particles, which should sup¬ 
press numerical heating. 

In Fig. 1^ we show the effect on the cooled gas fraction and 
on the number of formed galaxies with the asterisks. Quite interest¬ 
ingly, the resulting cooled fraction is very similar to that obtained 
for the MR and HR runs of the same cluster. As for the number 
of resolved galaxies, it is close to that obtained for the correspond¬ 
ing HR run, which has anyway a better gas-mass resolution. This 
increase in the number of identified galaxies inside the cluster can 
have two possible explanations. It is either due to an intrinsic in- 
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crease of the number of resolved halos where gas cooling takes 
place, or to an improved capability of small halos to survive dis¬ 
ruption in the cluster tidal field. To check this, we compare the 
values of ngai(< Rvii) inside the cluster with the corresponding 
number of galaxies in the region i?vir < R < 4f?vir around the 
cluster. Within Rvii, we identify 56 galaxies for the HR run and 53 
galaxies for the lo8 run. However, the difference between the two 
runs increases significantly in the outer cluster regions, with 123 
galaxies found for the HR run and 75 for the lo8 run. This shows 
that the better gas-mass resolution of the HR run in fact produces a 
larger number of galaxies in the field, while the improved accuracy 
in the gravitational force computation for the lo8 run compensates 
for this effect within clusters, where it plays a significant role in 
preventing tidal disruption of small galaxies. 

With respect to the cosmic star-formation history, the lo8 run 
(heavy solid curve in Fig.|3j shows a behaviour which is not strictly 
intermediate between MR and HR runs. The onset of star formation 
takes place at a higher redshift than in the HR runs, and then pro¬ 
ceeds in a more efficient way. Quite remarkably, the peak of star 
formation is even higher than for the highest resolution (VR) run. 
Finally, results of the lo8 run for the gas profiles are shown with the 
long-short-dashed curve in the bottom left panel of Fig.|^ Much 
like for the results on the cooled gas fraction, the resulting profiles 
are very similar to those of the MR and HR runs. 

Overall, the test on the degree of spurious gas heating from 
two-body relaxation confirms that a DM particle mass of at most 
~ 10® is sufficient to provide a reliable description of gas 

cooling within an already formed cluster of galaxies at low redshift. 
However, increasing the DM mass-resolution in order to reduce 
numerical heating of the gas from two-body encounters has a non- 
negligible effect on the number of resolved galaxies and on the star 
formation history, as a result of the better resolution of early struc¬ 
ture formation. The comparatively low computational cost of the 
gravity part in a hydrodynamical simulation with radiative cooling, 
and the availability of large amounts of memory in modem super¬ 
computers, may make it attractive to adopt a larger number of DM 
particles than gas particles in future simulation work. 

4.3 Grid vs. glass initial conditions 

As we discussed in Section|3 the initial Lagrangian particle distri¬ 
butions of the simulations of Set #1 have been realized as a grid, 
while those of Set #2 have been a glass. The gla ss-based technique 
to generate initial conditions (ICs; IWhitelHo^ aims at suppress¬ 
ing effects due to the regularity of the initial grid, which amplifies 
structure at the scale of the mean-interparticle separation. In the 
glass scheme, particle positions are initially generated randomly 
in the simulation box, but are then evolved backwards in time un¬ 
til they reach an amorphous, minimum energy configuration where 
each particle experiences only vanishingly small forces. The re¬ 
sulting irregular particle distribution lacks prefered directions and 
should be less affected by the symmetries that occur in the grid 
method. 

In order to check the effect of using either one of the two tech¬ 
niques to generate ICs, we have rerun the CL5 cluster with strong 
winds (SW), but this time starting from grid displacements. The 
results of this test are shown in Figure fTol where we plot the corre¬ 
sponding star-formation histories (left panel) and the radial profiles 
of gas properties (right panels). This comparison demonstrates that, 
at least at the resolution relevant for our simulations, the difference 
between using grid or glass ICs is very small. Looking at the de¬ 
tails of the comparison, it turns out that the grid-based run has a 


slightly higher star-formation at 6. This may be due to a con¬ 
tribution of the small-scale fluctuation modes around the Nyquist 
frequency, which should collapse more efficiently in the grid case 
(and part of this may be artificial), favoring somewhat earlier cool¬ 
ing at high-z. However, the star fractions within the cluster virial 
radius at « = 0 are /* = 0.17 for both runs. The number of identi¬ 
fied galaxies within the same region is 412 and 380 for the grid and 
glass runs, respectively. This appears to confirm that grid-based ICs 
show slightly higher power on small scales, which, in turn, gener¬ 
ates a slightly larger number of galaxies. As for the profiles of gas 
related quantities, we note that they also overlap quite closely. The 
only noteworthy difference is that the positions of merging sub¬ 
structures vary. However, such differences in orbital timing are ex¬ 
pected and common when different methods for the generation of 
ICs are used. 

The general result of this comparison is that, at least at the 
resolution relevant for our cluster simulations, the effect of using 
either grid or glass ICs is very small, in any case negligible with re¬ 
spect to other numerical effects that we have explored in this paper. 


5 CONCLUSIONS 

In this paper we have presented results from a large set of hy¬ 
drodynamical simulations of galaxy clusters, carried out with the 
Tree-l-SPH code GADGET-2. Our simulations include radiative cool¬ 
ing, star formation and energy feedback by a phenomenological 
model for galactic winds. The main target of our analysis has been 
the study of the stability of simulation results with respect to nu¬ 
merical parameters, such as mass resolution or gravitatiobal soft¬ 
ening length. We also considered different sources of numerical 
heating, and their interplay with the complex physical effects in¬ 
cluded. As such, our anal ysis also represents a validation study of 
our previous results te.g.. lBorgani et ai]l2004h . which were based 
on a large cosmological box simulated a relatively low resolution. 

Our simulated clusters span more than one order of magni¬ 
tude in collapsed mass and several decades in mass resolution. At 
the highest resolution, the mass of the gas particles is mgas — 
1.5 X 10^/i“^Mq, which allows us to resolve the virial region 
of a Virgo-like cluster with more than 2 million gas particles and 
at least as many dark-matter (DM) particles. Our main results are 
concerned with the effects of resolution on the properties of the 
stellar populations and on the intra-cluster medium of the simu¬ 
lated galaxy clusters. They can be summarized as follows. 

(a) In the absence of an efficient energy feedback, the fraction of 
cooled baryons steadely increases with resolution, reaching ~ 35 
per cent at the highest achieved resolution (VR runs in Table |3, 
with no indication for convergence. 

(b) Including feedback from galactic winds has the effect of stabi¬ 
lizing the stellar fraction inside clusters. Assuming a high efficiency 
for the supernova driven winds of order unity, we find that the frac¬ 
tion of cooled baryons converges already at a relatively modest res¬ 
olution, with an indication to decrease very slightly at the high¬ 
est resolution. This arises as a consequence of the self-regulation 
property of star formation and feedback. While improving the res¬ 
olution increases the star formation efficiency at very high redshift, 
this at the same time also provides a significant contribution to gas 
pre-heating which reduces star formation later on. The fraction of 
cooled baryons within R^ii lies in the range 12-18 per cent, with a 
decreasing trend with cluster mass. These values increase by about 
15 per cent when the normalization of the power-spectrum is raised 
from (78 = 0.8 to 0.9. 
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Figure 10. Comparison of the results for the CL5 cluster using both glass (solid lines) and grid (dashed lines) initial conditions. The left panel is for the redshift 
dependence of the star formation rate, while the right panels show the profiles of gas density, temperature, and entropy. Strong winds (SW) have been included 
in these simulations. 


(c) The feedback provides the necessary continued energy supply 
to keep gas particles of comparately low entropy and short cooling 
times in the hot phase of clusters, while without feedback these 
particles would cool down and drop out of the cluster atmosphere. 
As a result, the central gas density is higher in runs with feedback 
than in the runs with no galactic winds. The temperature profiles 
are shallower in the central cluster regions, while isentropic cores 
are much less pronounced, thus alleviating the discrepancy with the 
observed properties of the intra-cluster medium. 

A further series of tests presented in this paper concerns the 
effect of numerical heating. The main results from these tests can 
be summarized as follows. 

(a) Our Plummer-equivalent force softening of epi ~ 10 /i“^kpc 

at a mass resolution of moM — 5 x 10® scaled to other 

particle masses as epi oc represents a reasonable compro¬ 

mise between the need of preventing a spurious numerical heating 
of the gas, and the desire to resolve galaxies in the largest possible 
number of small DM halos. 

(b) Increasing the mass resolution in the DM component while 
keeping the gas mass resolution fixed, helps to reduce numerica l 
gas heating by two-body encounters jSteinmetz & Whitellio^ . 
We find that decreasing rriDM has a negligible effect on the frac¬ 
tion of cooled gas at 2 = 0 and on the radial profiles of gas proper¬ 
ties. However, it has a non-negligible effect on the star formation 
history, which occurs with enhanced efficiency at 2^ 3. This is the 
result both of a non-negligible numerical heating in small halos 
when gas and DM particles have rather different masses, and also 
of a better resolved dark matter content in the first generation of 
halos. 

In sum, our analysis helps to establish and delineate the regime 
of numerical reliability of the present generation of hydrodynam- 
ical SPH simulations of galaxy clusters. Even though quite tech¬ 
nical in nature, it is clear that systematic tests like the ones dis¬ 
cussed here are required to elucidate the often delicate interplay 


betweeen numerical effects and the physical model that is studied. 
Only when these effects are understood, the predictive power of nu¬ 
merical experiments can be fully exploited. At the same time, this 
understanding is also required to successfully push to new genera¬ 
tions of simulations with yet higher resolution. A reassuring aspect 
of the results discussed here is that even simulation models with 
complex and highly non-linear physical models for star formation 
and feedback can produce surprisingly robust results. This clearly 
is an encouragement for attempts to improve the fidelity with which 
the physics is represented in future simulation work. 
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